
#include "GeomDataLagrange.h"
#include "MyTime.h"
#include "BasisFunctionsLagrange.h"
#include "QuadratureUtil.h"
#include "ElementBase.h"

extern MyTime myTime;


GeomDataLagrange::GeomDataLagrange()
{
  bodyForceTimeFunction = -1;

  bodyForce.setZero();
}



void GeomDataLagrange::setNodalPositions(vector<myPoint>&  datatemp)
{
    nNode = datatemp.size();

    NodePosOrig.resize(nNode);
    NodePosCur.resize(nNode);
    NodePosNew.resize(nNode);

    specValCur.resize(nNode);
    specValNew.resize(nNode);

    myPoint  pt;
    for(int ii=0;ii<nNode;ii++)
    {
        pt = datatemp[ii];

        NodePosOrig[ii] = pt;
        NodePosCur[ii]  = pt;
        NodePosNew[ii]  = pt;
    }

    return;
}

void GeomDataLagrange::setNodalPositionsLocal(vector<myPoint>&  datatemp)
{
    int nNodelocal = datatemp.size();
    //
    NodePosOrigLocal.resize(nNodelocal);
    NodePosCur.resize(nNodelocal);
    NodePosNew.resize(nNodelocal);

    specValCur.resize(nNodelocal);
    specValNew.resize(nNodelocal);

    myPoint  pt;
    for(int ii=0;ii<nNodelocal;ii++)
    {
        pt = datatemp[ii];

        NodePosOrigLocal[ii] = pt;
        NodePosCur[ii]  = pt;
        NodePosNew[ii]  = pt;
    }

    return;
}



void GeomDataLagrange::updateNodesAfterPartition()
{
    myPoint  pt;
    for(int ii=0;ii<nNode;ii++)
    {
        pt = NodePosCur[node_map_get_old[ii]];

        NodePosOrig[ii] = pt;
    }

    for(int ii=0;ii<nNode;ii++)
    {
        pt = NodePosOrig[ii];

        NodePosCur[ii]  = pt;
        NodePosNew[ii]  = pt;
    }

    return;
}




void GeomDataLagrange::computeBasisFunctions1D(double uu, double *N, double *dN_dx)
{
    int  ROWS = 2, ii=0;
/*
    double** ders = new double*[ROWS], du_dx;

    for(ii=0;ii<ROWS;ii++)
       ders[ii] = new double[degree[0]+1];

    //HB_DersBasisFuns(degree[0], start, incr, uu, 1, ders);

    //  dx_du = J;
    //  du_dx = 1.0/dx_du

    //double  du_dx = 1.0/Jacobian[0];
    
    for(ii=0; ii<=degree[0]; ii++)
    {
        N[ii]     =  ders[0][ii];
        dN_dx[ii] =  ders[1][ii] * du_dx;
    }

    for(ii=0;ii<ROWS;ii++)
      delete [] ders[ii];

    delete [] ders;
*/
    return;
}



void  GeomDataLagrange::computeBasisFunctions2D(bool flag, int type, double* param, vector<int>& nodeNums, double *N, double *dN_dx, double *dN_dy, double& Jac)
{
    int  npElem = nodeNums.size();

    double  dN_du1[npElem], dN_du2[npElem];

    if(type == ELEM_SHAPE_TRIA)
    {
      LagrangeBasisFunsTria(npElem, param[0], param[1], N, dN_du1, dN_du2);
    }
    else  if(type == ELEM_SHAPE_QUAD)
    {
      LagrangeBasisFunsQuad(npElem, param[0], param[1], N, dN_du1, dN_du2);
    }
    else
    {
       //
    }

    // Gradient of mapping from parameter space to physical space
    double  xx, yy;
    MatrixXd  B(2,2), Binv(2,2);
    B.setZero();

    if(flag)
    {
       for(int ii=0; ii<npElem; ii++)
       {
          xx = NodePosCur[nodeNums[ii]][0];
          yy = NodePosCur[nodeNums[ii]][1];

          B(0,0) +=  (xx * dN_du1[ii]) ;
          B(1,0) +=  (xx * dN_du2[ii]) ;
          B(0,1) +=  (yy * dN_du1[ii]) ;
          B(1,1) +=  (yy * dN_du2[ii]) ;
       }
    }
    else
    {
       for(int ii=0; ii<npElem; ii++)
       {
          xx = NodePosOrig[nodeNums[ii]][0];
          yy = NodePosOrig[nodeNums[ii]][1];

          B(0,0) +=  (xx * dN_du1[ii]) ;
          B(1,0) +=  (xx * dN_du2[ii]) ;
          B(0,1) +=  (yy * dN_du1[ii]) ;
          B(1,1) +=  (yy * dN_du2[ii]) ;
       }
    }

    Jac  = B.determinant();
    Binv = B.inverse();

    // Compute derivatives of basis functions w.r.t physical coordinates
    for(int ii=0; ii<npElem; ii++)
    {
      dN_dx[ii] = dN_du1[ii] * Binv(0,0) + dN_du2[ii] * Binv(0,1);
      dN_dy[ii] = dN_du1[ii] * Binv(1,0) + dN_du2[ii] * Binv(1,1);
    }

    return;
}





void  GeomDataLagrange::computeBasisFunctions3D(bool flag, int type, double* param, vector<int>& nodeNums, double *N, double *dN_dx, double *dN_dy, double *dN_dz, double& Jac)
{
    int  npElem = nodeNums.size();

    double  dN_du1[npElem], dN_du2[npElem], dN_du3[npElem];

    if(type == ELEM_SHAPE_TETRA)
    {
      LagrangeBasisFunsTetra(npElem, param[0], param[1], param[2], N, dN_du1, dN_du2, dN_du3);
    }
    else if(type == ELEM_SHAPE_PYRAMID)
    {
      LagrangeBasisFunsPyramid(npElem, param[0], param[1], param[2], &N[0], &dN_du1[0], &dN_du2[0], &dN_du3[0]);
    }
    else if(type == ELEM_SHAPE_HEXA)
    {
      LagrangeBasisFunsHexa(npElem, param[0], param[1], param[2], &N[0], &dN_du1[0], &dN_du2[0], &dN_du3[0]);
    }
    else
    {
       //
    }


    // Gradient of mapping from parameter space to physical space
    double  xx, yy, zz;
    MatrixXd  B(3,3), Binv(3,3);
    B.setZero();

    if(flag)
    {
       for(int ii=0; ii<npElem; ii++)
       {
          xx = NodePosCur[nodeNums[ii]][0];
          yy = NodePosCur[nodeNums[ii]][1];
          zz = NodePosCur[nodeNums[ii]][2];

          B(0,0) +=  (xx * dN_du1[ii]) ;
          B(1,0) +=  (xx * dN_du2[ii]) ;
          B(2,0) +=  (xx * dN_du3[ii]) ;

          B(0,1) +=  (yy * dN_du1[ii]) ;
          B(1,1) +=  (yy * dN_du2[ii]) ;
          B(2,1) +=  (yy * dN_du3[ii]) ;

          B(0,2) +=  (zz * dN_du1[ii]) ;
          B(1,2) +=  (zz * dN_du2[ii]) ;
          B(2,2) +=  (zz * dN_du3[ii]) ;
      }
    }
    else
    {
      for(int ii=0; ii<npElem; ii++)
      {
          xx = NodePosOrig[nodeNums[ii]][0];
          yy = NodePosOrig[nodeNums[ii]][1];
          zz = NodePosOrig[nodeNums[ii]][2];

          B(0,0) +=  (xx * dN_du1[ii]) ;
          B(1,0) +=  (xx * dN_du2[ii]) ;
          B(2,0) +=  (xx * dN_du3[ii]) ;

          B(0,1) +=  (yy * dN_du1[ii]) ;
          B(1,1) +=  (yy * dN_du2[ii]) ;
          B(2,1) +=  (yy * dN_du3[ii]) ;

          B(0,2) +=  (zz * dN_du1[ii]) ;
          B(1,2) +=  (zz * dN_du2[ii]) ;
          B(2,2) +=  (zz * dN_du3[ii]) ;
      }
    }

    Jac  = B.determinant();
    Binv = B.inverse();

    // Compute derivatives of basis functions w.r.t physical coordinates
    for(int ii=0; ii<npElem; ii++)
    {
      dN_dx[ii] = dN_du1[ii] * Binv(0,0) + dN_du2[ii] * Binv(0,1) + dN_du3[ii] * Binv(0,2);
      dN_dy[ii] = dN_du1[ii] * Binv(1,0) + dN_du2[ii] * Binv(1,1) + dN_du3[ii] * Binv(1,2);
      dN_dz[ii] = dN_du1[ii] * Binv(2,0) + dN_du2[ii] * Binv(2,1) + dN_du3[ii] * Binv(2,2);
    }

    return;
}






